/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm.hxx"

#define UHM_ERROR_TOL 1.0e-5

int main (int argc, char **argv)
{
  FLA_Init();
  
  linal::Flat_ A1, A1_ext, A2, A2_ext;
  linal::Flat_ F1, F1_ext, F2, F2_ext;

  int
    verbose = 0,
    datatype = FLA_DOUBLE, 
    m = 3000, 
    n = 3000,
    offm_s = 1, 
    offn_s = 2, 
    offm_t = 3, 
    offn_t = 4,
    mm = 2, 
    nn = 2, is_erase = 0;

  A1.create( datatype, m, n );
  F1.create( datatype, m, n );
  A2.create( datatype, m, n );
  F2.create( datatype, m, n );

  FLA_Random_matrix( ~A1 );
  FLA_Copy( ~A1, ~F1 );

  
  if (verbose) {
    A1.disp("-A1-");
    F1.disp("-F1-");
  }
  
  A1.extract( A1_ext, mm, nn, offm_s, offn_s );
  A2.extract( A2_ext, mm, nn, offm_t, offn_t );
  
  if (verbose) {
    A1_ext.disp("- A1 EXT -");
    A2_ext.disp("- A2 EXT -");
  }

  {
    double t_base, t;

    int          m_B, n_B;
    int          rs_A, cs_A;
    int          rs_B, cs_B;
    char         blis_trans;

    rs_A     = FLA_Obj_row_stride( ~A1_ext );
    cs_A     = FLA_Obj_col_stride( ~A1_ext );

    m_B      = FLA_Obj_length( ~A2_ext );
    n_B      = FLA_Obj_width( ~A2_ext );
    rs_B     = FLA_Obj_row_stride( ~A2_ext );
    cs_B     = FLA_Obj_col_stride( ~A2_ext );

    double *buff_alpha = ( double * ) FLA_DOUBLE_PTR( FLA_ONE );
    double *buff_A     = ( double * ) FLA_DOUBLE_PTR( ~A1_ext );
    double *buff_B     = ( double * ) FLA_DOUBLE_PTR( ~A2_ext );

    FLA_Param_map_flame_to_blis_trans( FLA_NO_TRANSPOSE, &blis_trans );
    t_base = FLA_Clock();

    for (int i=0;i<1000;++i) {
      bli_daxpymt( blis_trans,
                   m_B, n_B,
                   buff_alpha,
                   buff_A, rs_A, cs_A,
                   buff_B, rs_B, cs_B );

      //FLA_Axpy( FLA_ONE, ~A1_ext, ~A2_ext ); 
      if (is_erase)
        FLA_Set( FLA_ZERO, ~A1_ext );
    }
    t = FLA_Clock() - t_base;
    printf("FLA_Axpy %lf\n", t);
  }
  if (verbose) {
    A1_ext.disp("- A1 EXT Axpy -");
    A2_ext.disp("- A2 EXT Axpy -");
  }

  F1.extract( F1_ext, mm, nn, offm_s, offn_s );
  F2.extract( F2_ext, mm, nn, offm_t, offn_t );

  if (verbose) {
    F1_ext.disp("- F1 EXT -");
    F2_ext.disp("- F2 EXT -");
  }
                               
  {

    double *buf_f1 = F1.get_buffer(), *buf_f2 = F2.get_buffer();
    double t_base, t;
    t_base = FLA_Clock();
    for (int i=0;i<1000;++i)
      uhm::merge(datatype,
                 buf_f1, F1_ext.get_offm(), F1_ext.get_offn(), m,
                 buf_f2, F2_ext.get_offm(), F2_ext.get_offn(), m,
                 mm, nn,
                 is_erase);
    t = FLA_Clock() - t_base;
    printf("Merging  %lf\n", t);
  }

  if (verbose) {
    F1_ext.disp("- F1 EXT Merge -");
    F2_ext.disp("- F2 EXT Merge -");
  }

  FLA_Axpy( FLA_MINUS_ONE, ~A2_ext, ~F2_ext );

  if (verbose)
    F2_ext.disp("- Diff -");
  
  linal::Flat_ norm_F1_ext, norm_F2_ext;
  norm_F1_ext.create( FLA_DOUBLE, 1, 1 );
  norm_F2_ext.create( FLA_DOUBLE, 1, 1 );
  
  FLA_Norm1( ~F1_ext, ~norm_F1_ext );
  FLA_Norm1( ~F2_ext, ~norm_F2_ext );
  
  printf("norm %lf, %lf\n", norm_F1_ext(0,0), norm_F2_ext(0,0) );


  norm_F1_ext.free();
  norm_F2_ext.free();
  
  A1.free(); A2.free();
  F1.free(); F2.free();

  FLA_Finalize();
  return 0;
}


